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ABSTRACT 


A  number  of  asymptotically  optimum  multidimensional  filtering  methods  are 
investigated  with  the  purpose  of  determining  filtering  techniques  which  require  rela¬ 
tively  little  computing  time  to  implement  with  a  digital  computer.  In  particular,  the 
asymptotic  properties  of  the  maximum -likelihood  and  minimum -variance  unbiased 
multidimensional  filters  are  investigated  in  the  sampled-data  case.  These  two  multi¬ 
dimensional  filters  are  shown  to  be  identical  since  they  are  both  based  on  a  conditional 
expectation.  In  addition,  the  martingale  property  of  conditional  expectation  assures 
that  the  asymptotic  properties  of  these  multidimensional  filters  are  well  defined. 

An  asymptotically  optimum  frequency  domain  synthesis  procedure  is  given  for 
two-sided  multidimensional  filters.  This  procedure  is  well  suited  to  machine  compu¬ 
tation  and  has  the  advantage  with  respect  to  the  exact  recursive  synthesis  method  of 
requiring  much  less  computation  time.  A  synthesis  procedure  for  physically  realizable 
multidimensional  filters  is  presented  which  is  based  on  a  factorization  of  rational 
spectral  matrices.  This  method  is  not,  however,  well  suited  to  machine  computation. 
An  interpretation  of  optimum  multidimensional  filtering  in  terms  of  frequency  wave- 
number  space  is  also  given. 
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I. 


INTRODUCTION 


It  has  been  recognized  that  a  considerable  reduction  of  seismic  noise  is  possible 
by  employing  multidimensional  filtering  for  seismic  arrays.1  Some  of  the  more  impor¬ 
tant  approaches  to  seismic  array  processing  are  the  maximum -likelihood  method,  the 
minimum -variance  unbiased  estimator  technique  and  multichannel  Wiener  filtering.  In 
general,  these  multidimensional  filtering  methods  form  a  single  output  waveform  which 
serves  as  an  estimator  of  the  unknown  signal  wi.ick  comes  from  a  fixed  direction. 

The  basic  assumption  in  the  analysis  of  multidimensional  filters  is  that  the  out- 
th 

put,  X^(t),  of  the  kl  seismometer  may  be  written  as 

Xk(t)  =  S(t)  +  Nk(t)  (1) 

where  S(t)  is  the  signal  waveform  which  is  assumed  to  be  the  same  in  each  seismometer 
and  Nk(t)  is  the  noise  present  in  seismometer  k,  k  =  1, . . . ,  K.  In  writing  Eq.  (1)  it  is 
assumed  that  the  azimuth  and  horizontal  velocity  of  the  event,  or  signal,  have  already 
been  determined  with  sufficient  accuracy  to  allow  the  signal  waveforms  from  each  seis¬ 
mometer  to  be  shifted  to  bring  them  into  time  coincidence.  In  most  applications,  the 
outputs  of  the  seismometers  are  given  in  sampled  form  in  which  case  Eq.  (1)  becomes 


.  +  N 
J 


k  =  1,...,K 
j  =  0,il,i2,  .  .  . 


(2) 


Only  the  sampled -data  multidimensional  filtering  problem  for  seismic  arrays  will  be 
considered. 
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It  can  be  shown  that  the  maximum -likelihood  estimator  is  the  same  as  the 
minimum  variance  unbiased  estimate  of  the  signal  when  the  noise  has  a  multidimen¬ 
sional  Gaussian  distribution.  In  addition,  the  multichannel  Wiener  filter  is  related 
very  simply  to  the  minimum  variance  unbiased  estimator.  The  synthesis  of  these  op¬ 
timum  multidimensional  filters  is  achieved  by  means  of  a  recursive  matrix  inversion 
algorithm  which  is  well -suited  to  computer  application.  The  purpose  of  the  present 
work  is  to  point  out  that,  asymptotically,  as  the  memory  of  the  filter  gets  large, 
there  are  synthesis  procedures  which  can  be  used  that  in  certain  cases  require  far 
fewer  computations  than  the  recursive  algorithm.  If  a  large  array  of  say  525  sensors 
is  to  be  processed,  then  these  approximate  synthesis  techniques  offer  a  considerable 
saving  in  computer  time  required  to  process  an  event  relative  to  the  exact  recursive 
procedure.  This  saving  will  be  shown  to  apply  only  when  the  filter  is  two-sided,  i.e. , 
non -physically  realizable.  This  restriction  poses  no  problem  in  the  application  of  the 
results  since  all  waveforms  are  recorded  on  magnetic  tape  and  non-physically  realiz¬ 
able  filters  are  readily  implemented. 

If  the  filters  are  specified  to  be  physically  realizable,  i.e. ,  one-sided  with  no 
output  before  any  input  is  applied,  then  the  synthesis  problem  becomes  more  compli¬ 
cated  than  in  the  two-sided  case  discussed  previously.  It  is  shown  that  a  spectral 
matrix  factorization  procedure  is  required  to  synthesize  the  filter  in  the  physically 
realizable  case  and  a  method  for  achieving  this  for  rational  spectral  matrices  is 
discussed.  Another  method,  due  to  Wiener  and  Masani,  is  also  described  which  is 
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valid  for  genera!  spectral  matrices.  Unfortunately,  these  spectral  matrix  factorization 
techniques  are  not  well  adapted  to  machine  computation. 

An  interpretation  of  optimum  multidimensional  filtering  in  terms  of  frequency 
wave-number  space  is  given.  The  structure  of  the  filter  in  frequency  wave -number 
space  is  also  presented. 


II.  DERIVATION  OF  THE  MAXIMUM -LIKELIHOOD  AND  MINIMUM -VARIANCE 
UNBIASED  ESTIMATORS 

The  derivation  of  the  maximum -likelihood  estimator  of  the  signal  requires  the 
assumption  that  the  noise  components  have  a  multidimensional  Gaussian  distribution. 
We  assume,  for  simplicity,  that  the  noise  components  have  zero  mean,  and  that  the 
covariance  matrix  is 


pkk/j,jP  "  E^NkjNk1j/ 


1  £  k,kj  £  K 
-  j,  ji  s  v  , 


(3) 


where  E  denotes  expectation,  and  it  is  assumed  that  the  estimator  is  to  use  2v  +  1 
samples  extending  in  time  from  -v  to  v.  Thus,  the  likelihood  function  can  be  written 


L  = 


-  |  (2vfl)  -  | 

(2tt)  |p  |  exp  { 


K 


4  2  2 
k,ki=i  j,jr'v 


(4) 
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where  jpj  denotes  the  determinant  of  the  matrix  p  which  is  a  matrix  of  K  x  K  submatrices, 
th 

the  j j ^  submatrix  has  the  elements  (j,jj),  k,kj  =  j,j1  =  with 

a  corresponding  notation  for  the  inverse  matrix  p  *  whose  elements  are  p.  /  (j»j,). 

kkj  I 

We  assume  throughout  that  the  matrix  p  is  positive  definite.  Differentiating  the 
logarithm  of  the  likelihood  function  with  respect  to  Sj^  and  equating  the  result  to  zero, 
we  obtain 


2 

k,kj=l  j=-v  1  J  J 


=  0, 


j  2  =  -v,  •  •  •  t  v  , 


(5) 


where  S  .  is  the  maximum -likelihood  estimator  for  S4. 
vj  J 

We  can  rewrite  Eq.  (5)  a  %  follows 


.  £  2  2  xkj  * 
js“V  1  k.k^l  1  k,kx=l  j=-v  1 


(6) 


J I  =  ”V» . .  • ,  v  • 

Let  us  define  the  (2v+l)  x  (2v+ 1)  matrix  G( j, j^)  in  terms  of  its  inverse  G  (j,  j^) 
whose  elements  are  defined  as 


-1/.  .  V 

a  (Jjp 


K 

Pui,  (j»ji)  » 

k,k=l  ttl  1 


j»j1  = 


(7) 


We  can  solve  for  S  .,  as 
vj 


V  -l  2  Mib')  ^  . 

J=-V  k=l 


(8) 


where 
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(9) 


0k(j|j,) 


K  v 

=  V  V 

Lj  Lj 

v 1  >rv 


(j. JA)  aUj.j') 


k  =  1 . K 

j  —  ”  V,  •  •  •  i  v  • 


It  should  be  no..  ,“d  that  a  different  set  of  weights  is  obtained  from  Eq.  (9)  for  different 
values  of  v  and  also  for  different  values  of  j',  which  can  vary  from  -v  to  v.  The  fact 
that  8k(  j  |  j!)  depends  on  v  has  not  been  brought  out  by  the  notation.  It  is  easily  seen 
from  the  quadratic  nature  of  the  logarithm  of  the  likelihood  function  that  the  solution 
given  in  Eq.  (9)  is  unique. 

We  now  consider  the  minimum -variance  unbiased  estimator  of  S.lf  denoted  by 

A 

S'  .. ,  expressible  as 
vj 


with  the  constraints 


(10) 


K 

Jj  \(jlj)  =  ® ii *  »  j  =  ,  (11) 

k=l  K  1J 

where 


j  =  j’ 

otherwise  . 


(12) 


The  variance  of  S'  is 
vj 


E(S^.)2  *  Z  Z  <4  Ojf)  \0lf>  ^  O.JP  • 

j.j^-v  k.k^l  1  1 


(13) 


Using  the  calculus  of  variations,  we  obtain  that  the  minimum -varia  nce  uuoiased 
estimator  has  weights  which  satisfy  the  system  of  equations 
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k  =  1, . 
j  =  -v, 


.K 

•  »v  » 


(14) 


K  v 

v  y 


>  6k1(jl,j,)  pkk1(j,jl)  +  ^vj 

r1  jrv  1  1 


=  0, 


where  the  X  are  2v+l  Lagrangian  multipliers  chosen  to  satisfy  the  constraints  given 
in  Eq.  (11). 

If  we  define 


0*K+l<j  IN  -  Xvj 


(15) 


Pk.K+l^^P  "  PR+l.k^^P  6k,K+l^ 

k  =  1 . K+l 

=  . . . 


(16) 


then  EqSv  (11)  and  (14)  may  be  written  as  a  single  set  of  equations  as  follows 


K+l 

E 

V1 


l 


jrv 


ek/jl^  *  Pkk,^*jP 


^k,  K+l  6jj* 


k  =  1 . K+l 

j  =  "V»  •  •  •  •  v  • 


(17) 


This  set  of  equations  is  equivalent  to  the  set  given  previously,  cf.  Reference  E ,  pp.  21- 

21.  The  reason  for  writing  the  equations  in  the  above  form  is  that  we  now  have  a 
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Toeplitz  matrix  of  submatrices  and  a  recursive  procedure  can  be  used  in  the  solution 
of  filter  weighting  coefficients. 

We  can  write  the  system  of  Eqs.  (14)  as 


6 


e^(jli')  ■ 


K  V  _ 

.2  Pkk.(j,jl)  Xvj  ’ 
=1  1  =-v  1 


k  =  1 . k 

j  “  ”  V*  *  •  •  ,  V 


Using  Eq.  (11)  in  (18)  we  get 


■  Z  Z  pkk  o.j,)\ 

k  kA=l  j1=-v  KK1 


5..,, 

JJ 


j  =  -V . . 


and  using  Eq.  (7) 


V 


Jl=-v 


=  ~  a(j,j’), 


j  ~  "V,  •  •  •  ,  v . 


Therefore,  we  see  from  Eqs.  (18)  and  (20)  that 

K  v 

9J,< j  b’)  =  Z  Z  ekk  ti.jp  a(j1,j’), 

V1  h=-v  1 


k  = 

j  =  -v . . , 


which  is  the  same  as  Eq.  (9).  Hence  the  maximum -likeliliood  and  minimum -variance 
unbiased  estimators  are  identical.  It  is  easily  seen  from  the  quadratic  nature  of  the 
expression  in  Eq.  (13)  that  the  solution  given  in  Eq.  (21)  is  unique. 

Using  Eqs.  (11),  (13),  and  (14)  we  obtain 


K  v 


e{svj.}  -  2  2  yjlm-y 


k-1  j— v 

—  f]  X  .6..,  •  -X  .,  *  a(j’,j')  , 

.  VI  11  VJ  » J  »  * 

j*-v 


j  •  -v,...  ,v 
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In  addition,  wc  have 


K 


E{S  S  .„} 
vj  vj 


■  l 


v 

V 

u 


k=l  i= 


@k(j I j’)  fa(j, 


i" 

J  /  j 


j=-v 


=  a(j \f). 


=  -V, 


»v  • 


(23) 


It  should  be  noted  that 


S  E{S*  .,}  , 


V*  —  V  • 


(24) 


This  follows  from  the  fact  that  the  minimum  value  of  a  quadratic  form,  subject  to 
certain  constraints,  in  a  v*  -dimensional  space  cannot  be  increased  if  the  dimension  of 
the  space  is  increased  to  v  ==  v'* 


HI.  THEORETICAL  JUSTIFICATION  FOR  THE  ASYMPTOTIC  APPROACH 

We  now  wish  to  show  that  the  use  of  filters  based  on  only  a  part  of  either  the 
past,  future,  or  past  and  future,  can  be  used  to  approximate  the  performance  of  filters 
based  on,  respectively,  the  full  past,  the  full  future,  or  the  Ml  past  and  future.  In 
order  to  do  this  it  will  be  convenient  to  introduce  new  random  variables,  as  indicated 
by  Rozanov,** 

Y,*.  =  K*- 
k*j  k*j 

Y..  =  X,*. -X..  k^k*,  k=  1,...,K  (25) 

kj  ^  j  =  0,±l.±2,... 
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Let  us  denote  the  set  of  random  variables  Y  ,  k  j-  k*  j  =  -v . .  by  Y  .  It  should 

kj  v 

be  noted  that  the  random  variables  Y  ,  k  i  k*,  do  not  depend  on  S.,  i.e. , 

kj  j 


Ykj  -  Nk*rv 


k  i  k*, 


(26) 


and  that 


Y  =  S.  +  NU.  . 
k  J  J  k*j 


(27) 


It  therefore  follows  that  the  minimum -variance  unbiased  estimator  of  S.t  may  be  written 


as 


S  .,  =  X  #1,  N  i 
VJ  k*j  vk*j 


(28) 


where  N  ,  is  a  linear  combination  of  the  Y,  .’s,  k  ^  k*.  j  -  -v, . . . .  v,  i.e. ,  N  , 

vk*j  kj  vk*j 

A 

depends  only  on  Y^.  When  written  in  this  form,  it  follows  easily  that  is  given 

by 


N  =  E{N  ,  ,.,|Y  }, 

Vk*j’  \ac  j  v 


(29) 


i.e. ,  N^*j,  is  the  conditional  expectation  of  N^*.,  with  respect  to  Y^,  cf.  reference  6, 

A 

pp.  150-155.  However,  N  is  a  martingale  since,  cf.  reference  6,  pp.  91-94, 

VK  J 

E,Vl1'|Vv}  ■  E{E|V|Yv+ll|Yv} 


■  E<V|Yv}  ■  Vr  • 


<30) 
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■  I  *  — niwpiii  .nwnitti  niiiiiiniiiiT-  iinf 'in ii'Wim 


Thus,  it  follows  from  a  martingale  convergence  theorem,  cf.  reference  6,  p.  167, 
Theorem  7.4,  and  also  pp.  560-562, 


l.  i.  m. 
v  -* 00 


N 

vk*j*  * 


(31) 


for  fixed  j*.  If  we  wish  to  consider  sequences  of  physically  realizable  filters  for  which 
j'  =  v,  we  obtain  in  a  similar  manner 


l.  i.  m. 
v  00 


N 


k*v 


(32) 


and  for  filters  using  future  values,  j*  =  -v,  and 


t.  i.m. 
v-*  ® 


N 


vk*,  -v  ' 


(33) 


The  minimum -variance  unbiased  estimators  for  the  above  three  cases  are, 
respectively, 


Sj,  =  l i-  na.  Svj,  ^*j*  ^k*j'  *  (^ 


A 

s 

00 


l.  i.  m. 

v  ” 


A 


(35) 


A 


s 

-  00 


l.  i.  m. 

V-*  00 


S 

V,  "V 


Xk*,-»  Nk*, -00 


(36) 


Thus,  we  have  shown  that  filters  based  on  a  large  part  of  either  the  past,  future,  or  past 
and  future,  can  approximate  the  performance  of  filters  based  on,  respectively,  the  full 
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past,  the  full  future,  or  the  full  past  and  future.  It  is  easily  seen  that  all  results 

remain  valid  when  k*  is  any  integer  between  l  and  K. 

We  have  already  shown  that  the  maximum -likelihood  estimator  for  S.,  is 

J 

identical  with  the  minimum -variance  unbiased  estimator  for  S^,,  cf.  Eqs.  (9)  and  (21). 
However,  at  this  point  it  is  extremely  simple  to  show  that  this  is  true.  The  joint 
probability  density  of  X^*.,,  Y  ^  is 


<Vr  V  ■  psj.<WYv>p(V 


(37) 


since  is  independent  of  S^,.  In  addition,  we  have  for  the  Gaussian  multivariate  case 


PS 


.<Vj.lV  ■  W  1/2PtXk...|Yvjexp{— |  l  i2i 


(38) 


where  j(X 


vy 


Y^)  is  independent  of  S^,  and 


=  E{tVj' 


V 


1-S 


vj 


(39) 


The  maximum -likelihood  estimator  for  S^,  is  obtained  by  differentiating;  with  respect 
to  Sj,,  the  probability  density  function  in  Eq.  (37),  or  equivalently  that  in  Eq.  (38), 
from  which  we  get,  when  using  Eq.  (39), 


(40) 
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which  agrees  with  Eqs.  (28)  and  (29).  Thus,  all  asymptotic  results  derived  for  the 
mi.iimum -variance  unbiased  estimator  remain  true,  of  course,  for  the  maximum - 
likelihood  estimator. 


IV.  APPROXIMATE  FRECUENCY-DOMAIN  SYNTHESIS  PROCEDURE  FOR  TWO- 
SIDED  FILTERS 

An  approximate  frequency-domain  synthesis  procedure  will  now  be  presented 
for  maximum -likelihood  filters  which  use  a  large  part  of  the  past  and  future.  Such 
filters  which  use  the  past  as  well  as  future  values  will  be  termed  two-sided  filters. 
We  begin  by  assuming  the  noise  is  wide-sense  stationary  and  by  letting  v  =  ®,  j'  =  0, 
so  that  Eqs.  (13)  and  (14)  may  be  written  as 


E{S2} 

o 


limit  E{S  } 
vo 


rr 

|* 

J 


2  f  (x)  A*(x)  A  (x) 
j,k=l  JK  J  * 


dx 

2tt  * 


(41) 


and 


where 


;  j  ,  <x)  Ak<x)  €-,te  %  +  h 

■tt  k-i  J 


=  0, 


—  0,  ±1,  ±2, . 

j  =  1».  •  •  »K, 


(42) 


Ak(x)  =  £  0k(j|O)€1JX. 

j=-C0 


k  "*  1,*.., X, 


(43) 
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/#\  p  y  .  -i4x  dx 

=  J  Vx)  €  55  • 

—  TT 

and 

j,  k  =  1, . . . ,  K 

4—0,  il,  i  2, . . . 

V»  ■  Jj.  €‘"x  ■ 

j,  k=l, . . . ,  K  , 

(44) 

is  the  sampled  cross  power  spectral  density  function,  x  =  ujT,  T  is  the  sampling 

interval.  We  have  from  Eq.  (43)  that 

ek(j|o,  -  j\<x>  |  , 

-TT 

k  =  1,...  ,K 
j  =  0, il,  ±2, . . . 

(45) 

The  constraint  Eq.  (11)  becomes 


K 

E  V*> =  1  •  (46> 

k=i  K 


If  we  let 


AM  -  E  • 


(47) 


k=-°° 


then  Eq.  (42)  may  be  written  as 


TT 


J  € 


-i-tx 


-rr 


k=l 


f -k<x)  Ak(x>  +  A(x> 


dx  =  0  , 


4  =  0*  il,  ±2,  •  • 
j  =  1,...,K  . 


(48) 


According  to  Eq.  (48),  the  Fourier  coefficients  of  the  quantity  in  the  brackets 
in  E  )  must  all  be  zero.  It  follows  that  the  quantity  itself  must  be  zero  so  that 
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(49) 


K 

Z  fit(x)  A,<x)  +  A(x)  =  0  , 

k=l  J 


-tt  2  x  £  n 

j  ”  i « » j k  * 


Thus 


Ak(x) 


K 

Z  qki<x> 

JE±_  J 

•Zvx> 

J>k=l 


K 

A(x)  =j  ^  q(x) 

L  j » k= 1 


-1 


k  =  1, . . .  ,K  ,  (50) 


(51) 


where  {q^OOl  is  the  inverse  of  the  spectral  matrix  {f^(x)},  j,k  =  1,. . .  ,K.  We  note 
that 


Vx>  ■  qkVx)  • 


(52) 


since 

V(x>  =  fkj(x)  ■ 

and 

qjlr(x)  --  M*k(-x)  . 

since 


(53) 


(54) 


V(x>  =  ^('x> 


(55) 


Therefore 
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A  (x)  =  A*(-x)  , 
k  K 


k  =  1,.. .  ,K 


(56) 


We  obtain  from  Eqs.  (41)  and  (50) 


rr 


Ms2}  -  r  £ 

1  o  J  2n 

-  TT 


K 


■  k  i  V0 

j,k=l 


(57) 


The  filter  weighting  coefficients  are  obtained  from  Eq.  (45)  as 


rr 


Q.(jlO)  =  J  ~  [Re  A  (x)cos  jx  +  Im  A  (x)  sin  jx],  k=l,...,K 


—  TT 


2tt 


j  =  0,±1,±2,, 


(58) 


It  is  easily  seen  that  the  constraint  conditions  are  satisfied  since 

K  TT  , 

Z  ©k( J  |0)  =  J  cosjx  ^ 

k=l  -TT 

=  6.  •  (59) 

jo 

The  filter  weights  given  by  Eq.  (58)  would  have  to  be  obtained  by  inverting  the 
spectral  matrix  at  all  points  on  the  frequency  axis  between  -rr  and  n.  This  is  clearly 
impractical  so  that  an  approximation  would  have  to  be  used  in  practice.  Such  an 
approximation  is  most  easily  obtained  by  approximating  the  integral  in  Eq.  (58)  by  a 
finite  sum 
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{||1|jl;j|,|,|j,||||]|  j|j|jj||  jf|,|  j  j  i  j  I  i : !  ]  1  iijHiilHhilii  mi:  1 1  Ilf  .ulilllHlilllli  jl-llli*  ‘ 


(60) 


ek(j|0) 


-  —  y  Rc  A.  (I  -)  cos  }l-  +  Im  A  (-1  -)  sm  \l  -  , 

2v  #  u  , ,  k  v  v  k'  v  *  v 

i=-v+l 

k  =  1 . K 

j  =  -v . v 


The  symmetry  properties  of  A^OO  expressed  in  Eq.  (56)  enable  us  to  simplify  the  above 
equation 


v- 1 

9.  (j|0)  =  —  f  A.  (0)  +  (-1)J  A  (rr)  +  2  y  Re  A  At  -)  cos  j-l-  +  Im  A  (l -)  sin  j-t  -  ] 
k  2v  -  k  k'  '  ,H  k  v  J  v  k  v  J  v 


t=l 


k  1( .  i. . ,  K 
j  =  “V»  -  •  •  *  v 


(61) 


TTie  constraint  equations  are  still  satisfied  since 


K 


Z  vil«>  -  * 


k=l 


_L 

2v 


v-1 


(-l))  +  2(i+  Yj  COS  jt^> 
•fc-i 

(-1)*  +  Sln<v-I>^  v 


.  1  .  IT 

sui-J- 


-  [+1  +  2V-1]  =  1, 


j  =0 


2v 


(-DJ  - 


/If  •_  \ 

sm(jj--jTT) 

.  1  .  IT 

sm  -  j  - 
2  v 


=  h  K-i)1  -(-I)3]  =  0  , 


j  , 


(62) 


where  we  have  used  the  identity 
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n 

V 

u 

k=l 


cos  ku 


x  sin(n+-)u 


2  .  1 
Sln  2  U 


(63) 


It  is  easy  to  see  that  as  v  the  estimator  based  on  the  approximate  filter  weights 

A 

given  by  Eq.  (61)  converges  in  the  mean  to  S^. 

In  order  for  the  asymptotic  results  to  hold,  we  must  have 


2vTW  >  1 


(64) 


where  T  is  the  sampling  interval,  and  W  is  the  bandwidth  of  the  major  spectral  peak 
of  the  f.^(x)  and  is  assumed  to  be  roughly  the  same  for  all  j,k.  It  has  been  found  ex¬ 
perimentally  that  for  microseismic  signals  T  should  be  about  0. 1  seconds  so  that 
W  >  5/v  cps.  For  21  filter  points,  v  =  10,  so  that  W  >  ^-cps.  It  appears  from  pre¬ 
liminary  spectral  analysis  that  W  is  approximately  j  cps,  or  slightly  smaller,  so  that 
it  is  difficult  to  say  whether  the  time -bandwidth  product  given  in  Eq.  (64)  is  sufficiently 
large  for  the  asymptotic  results  to  hold  for  21  filter  points.  In  general,  the  rate  with 
which  the  asymptotic  results  are  achieved  would  have  to  be  determined  experimentally 
by  means  of  computer  calculations. 

We  now  compare  the  asymptotic  results  obtained  for  the  maximum -likelihood 
filter,  with  those  for  the  Wiener  filter.  The  Wiener  filter  functions  are  given  by  (see 
Reference  2,  p.  10), 


Hk(x)  =  Ak(x)  G(x)-A(x)  * 
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. 


•  lere  G(x)  is  the  assumed  spectral  density  function  for  the  signal,  when  the  signal  is 
ta,.  _n  to  be  a  stationary  random  process.  These  filters  are  the  seme  as  the  maximum  - 
likelihood  filters  multiplied  by  the  common  filter  function 

G(x) 

G(x)  -  A(x)  ' 

The  filter  weighting  coefficients  are 


\  =  f  G 


7  Go«e-U“  dx 


(x)  -  A(x)  2rr  * 


which  may  be  approximated  as 


i  #  Go<i  v> 

"k'Str.i  ;;.n.  cos)kv 


k  =  0,±l,  +  2, . . .  (65) 


k  ~v,  •  •  • ,  v  (66) 


We  now  compare  the  computing  times  required  b>  die  exact  recursive  and 
approximate  frequency  domain  procedures.  It  is  easily  seen  that  the  exact  recursive 
procedure  requires  approximately,  cf.  Reference  4,  pp.  98-105, 

10v^  kt  +  a)  seconds  ,  (67) 

where  the  amount  of  computing  time  required  to  invert  a  K  x  K  matrix  is  taken  to  be 

3 

K  (n  +  a),  m  and  a  are  the  multiplication  and  addition  times,  respectively,  in  seconds. 

If  v  and  K  are  large,  the  approximate  frequency  domain  synthesis  procedure  requires 
essentially  the  inversion  of  (v  +  1)  K  x  K  Hermitian  spectral  matrices  plus  a  Fourier 
transform  operation.  The  matrix  inversion  requires 
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(v  +  1)  *  4  *  4K^  (m  +  a)  —  2vK^  (p.  +  a)  seconds 


and  the  I  ourier  transform  operation  requires 


4v  K(u  +  a)  seconds 


for  a  total  of 


2vK^  (|j  +  a)  [  1  +  2  ~  ]  ^  2vK^  (|_t  +  a)  seconds  , 

T/* 


f»  i. 

2v 


In  Eq.  (68)  the  factor  4K  is  required  for  the  inversion  of  matrices  with  complex 
elements  and  the  factor  of  j  enters  due  to  the  symmetry  of  the  matrices.  Thus,  the 
approximate  frequency  domain  synthesis  procedure  requires  5v  times  less  computing 
time  than  the  recursive  procedure. 

A  convenient  way  to  estimate  the  f_k*s  is  to  transform  the  estimate  of  the 
correlation  coefficie.  ts,  i.e. , 


fjk« 


2v 

-  Z  o 

l=-2v 


-  ^ 


j,k  =  1, . . .  ,  K  , 


where  p  (l)  is  the  estimated  correlation  coefficient  obtained  in  the  following  way, 
J* 


i  2 

L  0  -  is  L 
O^i+l^L 


N.  . . f  N.  .  , 

j,i+t  k,i 


j,k  =  1,...,K 
K/  ~  0,1 . 2v  , 


where  L  is  the  number  of  samples  used  to  estimate  the  correlation  coefficients.  This 


method  for  estimating  f  (x)  leads  to  estimated  spectral  matrices  which  are  nonnegative- 

JK 

A  A 

definite  for  all  frequencies,  -n  <  x  ^  tt.  Note  that  p  (-1)  =  p  .(£),  l  =  0, 1,. . .  ,2v. 
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H'Wii»lii.;wl|ilMlt  IfflBMI!  WW#5 


The  estimated  covariance  .  atrix  is  nonnegative -definite  since 


K  v 


2  2  P.,(m,n)  a.  cl* 

. T  i  jk  lm  kn 

J,k=l  m, n=-v  J  J 


.  K  v 

■i  E  2 

j,k=l  m,n=-v 


V 


A  N,  .  N. .  a.  a,* 

0<i£L  J ,  i+m  -i:  ki  jm  Hen 


0  £  i+m-n  L 


,  K  v 

=  i  y  v 

t  Li  lj 


2  N.  N.  a.  a* 

“j,k=l  m,n=-v  OStan  Jll+m  k>H”  )m>n 
0  S  i4n  ^  L 

L  K  v 

=  r  2  2  2  e(i+m)  e(i+n)  N.  N  a.  a.* 

L  i  i  u  j.i+m  k,i+n  p  Ten 

1=0  j,k=l  m,n=-v  *  J 


f  t 

i=0 


K  v 


2  2  e(x + m) a-  N- . . 

u  im  1,1+m 

1=1  m=-v 


s  0  , 


where 


e(i  +  m)  =  1,  O^i  +  msL 


=  Of  otherwise  . 


a 

We  now  show  that  the  estimated  spectral  matrix  {f .^(x)}  is  nonnegative- 


definite  for  all  x,  by  noting 


V>- 11'-*’  voe 


i-Cx 


so  that 


JL-  l 

)+  1  U 


2v  +  1  u  Pjk(”.")« 
m,n=-v 


i(m-n)x 
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K 

V 

/  , 


a 


K 

V 


v 

V 


j,k=l 


.  a*  f..  (x)  =  )  a,  a,*  — - 7  ■ 

]  k  jk  .  r1  ,  j  Tc  2v  +  1  ^ 

J,k=l  n,m= 


■v 


n.k(m,n) 


ci(ni  -n)x 


i _  V  „  V  ,-‘(m-n)x  I  1  y 

>+1,4.  “A  ^  €  ,L  „4 


2v+1j.^rjntn,.^  ‘  lL  OSUL  Ni.^A 


0  -  i+m-n  —  L 


_ i _  y 

L(2v+1)  £Q 


2  2  ea+m>aeiraxN 

j=l  m=-v  J 


^  0  . 


Ordinarily,  the  amount  of  computer  time  required  to  perform  the  Fourier  transforma¬ 
tion  in  Eq.  (7 1)  would  be 


(2*4v)  (v+  1)  ~  (|a+a)  -  4v2  K2  (ia+a)  seconds 


and  this  time  should  be  added  to  that  in  Eq.  (70).  If  v  and  K  are  of  the  same  order  of 
magnitude  then  the  above  computer  time  is  comparable  to  that  in  Eq.  (70).  However, 
it  is  possible  to  compute  the  Fourier  transform  in  Eq.  (71)  in  less  than,  cf.  Reference 

7, 

K2 

(8v  log2  4v)  —  (m  +  a)  —  4v  K2  log2  4v  (m  +  a)  seconds. 

Thus,  if  v  is  large  so  that 


log9  4v 

— - -  «  1, 
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a  considerable  saving  in  the  amount  of  computer  time  required  to  perform  the  Fourier 
transform  in  Eq.  (71)  is  achieved.  In  addition,  if  we  have 

}  _ K _  >;> 

2  log0  4v  >:>  * 

Am 

then  the  amount  of  computer  time  required  for  the  Fourier  transformation  becomes 
negligible  compared  to  that  in  Eq.  (70),  and  may  therefore  be  neglected. 
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V.  SYNTHESIS  OF  PHYSICALLY  REALIZABLE  FILTER  BY  MEANS  OF  A  SPECTRAL 
MATRIX  FACTORIZATION  METHOD 


We  now  consider  the  synthesis  of  the  filter  for  the  physically  realizable  case, 
j'=  v,  v~* 00 .  In  this  case  we  have 


Ak(x) 


00 


2 

j«0 


ek(jb) 


and 


tt  K 

J*  Zl  {  E  f.k(x)Ak(x)  +  A(x)}  dx  =  0 

TT  lr—  1  * 


L  =  0, 1,  2,... 
j  =  1, . . . ,  K. 


(72) 


(73) 


Let  us  define  z  =  € 


ix 


and 


K 

*.(z)  =  2  A(z)  +  A(z)  ,  j  =  l,...,K.  (74) 

J  k=!  k 


We  have  from  Eq.  (72) 


A. (z)  =  E 
j=0 


0k(j|°°)  . 


k  =  1, . . . ,  K  .  (75) 


Since  A  (z)  is  a  power  series  in  ascending  powers  of  z,  A  (z)  must  be  analytic 

It  It 

in  the  unit  circle  of  the  z-plane.  In  order  to  have  tk(z)  satisfy  Eq.  (73)  we  must  have 


k=  1, 


(76) 
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so  that  |'^(z)  is  analytic  outside  and  on  the  unit  circle  of  the  z-plane.  If  we  subtract 
the  first  equation  in  (74)  from  all  the  others  and  substitute 

K 

A  (z)  =  1  -  £  A  (z)  ,  (77) 

k=2 

we  obtain  a  new  system  of  equations 

K 

Z  [f11(z)  +  fik<z)-fj,(z)-flk(z)l  Ak(z)  =  fu<z)-fjl(z)  +  i1(z)-}.(z»  , 
k=2  J  J  J  J 

j  =  2 . K  .  (78) 

The  functions  i|f  j(z)  -  i|l(z)  ,  j  =  2,. . .  ,K  are  analytic  outside  and  on  the  unit  circle. 

Thus,  we  have  a  system  of  (K-l)  equations  in  (K-1)  unknowns  from  which  we 
can  solve  for  A2(z), '  *  *  We  then  obtain  A^(z)  from  Eq.  (77).  The  system  of 

equations  in  (78)  may  be  written  in  matrix  notation  as 

A(z)  f(z)  =  h(z)  +  t|i(z).  (79) 

The  matrix  f(x)  is  easily  seen  to  be  a  spectral  matrix.  It  will  be  assumed  that  f(z)  has 
a  spectral  matrix  factorization 

f(z)  =  P(z)P’(z'1)  ,  (80) 

where  the  matrices  P(z),  (P(z)]  *  have  matrix  Laurent  expansions  on  jzj  =1  with  no 
negative  powers  of  z  and  F*(z)  denotes  transpose.  We  have 
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AfcJPfeJFfc-1)  =  h(z)  +  f(z). 


(81) 


and 

A(z)P(z)  =  h(z)  [P!(z_1) )  1  +  i(z)  [  P* (z ~ 1 )  ] ~ 1 


(82) 


The  matrix  on  the  left-hand  side  of  Eq.  (82)  has  a  matrix  Laurent  expansion  with  no 
negative  powers  of  z  which  is  assumed  to  converge  in  some  annulus  containing  the  unit 
circle.  The  second  term  Oi  he  right  has  only  negative  powers  of  z.  Equating  coefficients 

o 

we  obtain,  cf.  Whittle,  pp.  67,  100, 

A(z)  =  {h(z)  IP*(z-l)]_1}+  IP(z)]'1  ,  (83) 


where  the  operation  {  }+  indicates  that  only  the  non-negative  powers  of  z  in  the  Laurent 

expansion  of  the  matrix  within  the  braces  are  to  be  retained.  Equation  (83)  represents 
the  complete  solution  to  the  synthesis  problem  for  multidimensional  physically  realizable 
filters. 

It  is  now  necessary  to  show  how  the  factorization  in  Eq.  (80)  can  be  obtained. 

A  procedure  for  obtaining  the  spectral  matrix  factorization  for  rational  spectral 

g 

matrices  has  been  given  by  Whittle,  pp.  101-103,  and  is  similar  to  that  used  in  the 
one -dimensional  case,  cf.  also  Rozanov.^’ ^ 

As  an  example,  let  us  consider  that  K  =  2  and 
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f  (z)  +  f  (z)  —  f  (z)  —  f  (z)  *  — — — X 

*ir  j  22l  '  12*'  2V  '  (1-oz)  (1  -az'1) 


--if  .  isl  <  1  • 


<85) 


We  have 


FW  -  f'1Z-  , 

1  -  oz 


(86) 


and 


A2(z>  = 


1  T  .1  ir 

J1-0Z)  (l-pz-l)J  1  - 


-  az 
@z 


(87) 


Since  |a.j  <  1,  the  Laurent  expansion  of  l/(l-ctz)  may  be  made  in  positive  powers  of  z 

and  will  converge  in  the  circle  |z|  <  »  which  encloses  the  unit  circle.  The 

Laurent  expansion  of  l/(z-0)  must  be  made  in  negative  powers  of  z,  i.e. , 

"1  “2 

z  +  0z  +  ... ,  and  will  converge  in  the  annulus  1 0 1  <  |z  |  <  °°,  which  also  encloses 
the  unit  circle.  Th.  s,  in  order  to  perform  the  required  operation  on  the  term  in  the 
brackets  in  Eq.  (87)  we  perform  a  partial  fraction  expansion  and  neglect  the  term 
l/(z-0),  to  obtain 


A2(z)  =  l-a0  l-0z 


(88) 


a  -  JL_  z(og-i)  -  a 

Al{  '  1-00  1-Bz 


(89) 
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We  now  wish  to  give  an  example  i.i  which  we  will,  compare  the  performance  of 


the  physically  realizable  filter  with  that  of  the  two-sided  filter.  For  simplicity  we  con- 
sidei  that  K  =  2  and 

f  (z)  =  - - - z~r  ,  0  <  s.  <  1 

(l-az)(l-az  ) 

f9  (z)  = - - - 7-  ,  0 <  8<  1 

(i  -  ez)  (i  -  az  L) 


fJ2(z)  =  =  0  * 


so  that 


fll(z)_f21(z)  = 


1 


(1-oz)  (1 


az  *) 


fn(z)  +  f22(z)  -  f12(z)  -  f21(z) 


(1  -  qz)  (1  -az1)  +  (1  -  8z)  (1  -  9z  S 
(l-az)(l-az  S  (1  -  gz)  (1  -  Bz  J) 


In  addition,  let  us  define 

(1  -  8z)  (1  -  8z  A)  =  (i  -  az)  (1  -  az  A)  +  (l  -  0z)  (1  -  8z  *) 


so  that 


0  =  a  +  0 
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and 


0^  =  i  +  a2  +  er3  . 

Thus,  a  and  0  must  satisfy  the  equation 

2a  0  =  1 


and 


e 


Hence  a  and  0  must  lie  in  the  open  interval  ( -  ,  1),  0  is  in  the  open  interval 
<VT,  -),  and  we  have 


fll(z)+f22(z)_f12<z)"f21<z) 


ffp-e^z)  (l-e'V1) _ 

(1  -  az)  (1  -  az  2)  (1  -  0z)  (1  “  0Z  *) 


Therefore,  we  can  write 


P(z)  = 


9(1  -  01  z) 

(1  -cu.)  (1  -  0z) 


and 


A2(z) 


CD 

h- 

1 

CD 

1 

H— * 

N 

_ 1 

(1  -  az  l)  (1  -  0z  l) 

(1  -  az)  (1  -  0z) 

9  (1  -  az)  (1  -  az  ^{l-e  *z  *) 

a 

l  -  az  * 
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1  z  -  2a 


> 


A^z) 


1 

2 


-1 

z  -a 
z  -  G 


We  have  from  Eq.  (13) 


I  A.(X)A*  (x)fjk(x)  | 
J*  *-1 


S  ^  E  A  (z)  Ak(z-‘)f(z) 
unit  j,k=l  J 


dz 

z 


circle 


K 


=  2  Residues  \  £  -  A  (z)  A  (z  *)  f  (z) 


inside 

unit 

circle 

=  Res  < 
z=S_1 


j,k=l  2  j  k  jk 


z  -a 


-1 


+  ? 


z  -  2a 


4a(z  -  0)  (9z  -  1)  (1  -  az)  2  (z  -  0)  (0z  -  1)  (1  -  0z) 


1 

0s  -1 


[a2  + 


l 

4a2 


]  • 


For  the  two-sided  filter,  we  obtain 
qn(z)  =  (1-az)  (1  -  az  *) 
q22<z)  =  (*  “  Bz)  (*  “  Bz  S 
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qi2<z>  =  q2l^  =  0 


2  l"1 

E  qjkfe)  - - - - — 

j,k=l  J  (1  -  6z)  (1  -  0z  ) 


Using  Eq.  (57)  we  l.av  •, 


A0  rr  rK  i-l  , 

=  E  y»  | 

-n  ],k=l 


Thus 


*  2  V>  f 

unit  j,k=l 


circle 

=  I  Res  I  \  S  qik<z>  [ 

inside  |_j,k=l  J  | 

unit  *  ' 

circle 


z^q-1  (l-9z)(z-0) 

1 

0^-1 


E{S2| 

c 


2  ^  i 
a  4a2 


The  minimum  value  of  the  above  ratio  is  unity  and  occurs  when  a  =  1/a/2  .  The  maxi¬ 
mum  value,  in  the  permissible  range  for  a,  is  5/4  and  occurs  when  a-^ora"  1.  Thus, 
there  can  be  a  loss  in  noise  variance  reduction  of  between  zero  and  approximately  1  db 
by  using  the  physically  realizable  filter  rather  than  the  two-sided  filter. 
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Another  method  for  achieving  the  spectral  matrix  factorization  which  is  valid 
for  general  spectral  matrices  is  due  to  Wiener  and  Mason).  Suppose 


0<  ml  ^  1 1 f (x)  |1  s  MI<  ®  , 


where  the  norm  ||a  j|  of  the  matrix  a  =  {a  }  is  equal  to  ) 

J  J.k 


2 

iajk  1  .  and 


f(x)  =  I+g(x)  , 

llsCx)  ||<  J  • 


Under  these  conditions  we  have 

f(x)  =  b_1(x)  C1/2[b_1(x)  C1/2]'  , 

where 

b(x)  =  I-g  +  (g  g)  ~  [<g  g)  g]  +  ... 

OC 

the  matrix  g  (x)  is  obtained  from  the  matrix  g(x)  =  T)  g  S  by  omitting  the 

0  ..  k=-»' 

IX  -  n  1KX 

positive  powers  of  €  ,  g  (x)  =  Si.  €  .  and  the  matrix  C  is  defined  by 

k= 

C  =  b(x)  f(x)  b*(x) . 

-1  1/2 

Thus,  the  factorization  of  f(x)  is  achieved  in  terms  of  b  (x)  C  .  It  should  be 
noted  that  C  is  a  constant  nonnegarive -definite  matrix  which  is  independent  of  x  so  that 
its  square  root  may  be  obtained  in  the  usual  way,  i.  e. ,  by  diagonalizing  C  with  a 
unitary  matrix  and  taking  the  square  root  of  the  resultant  diagonal  matrix.  Unfortunately, 
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b(x)  can  be  obtained  only  as  an  infinite  sum  of  matrices,  the  rate  of  convergence  of 
which  is  difficult  to  determine. 

It  is  readily  appreciated  that  the  techniques  for  synthesizing  physically 
realizable  filters  by  means  of  a  spectral  matrix  factorization  method  for  rational 
matrices  are  impractical  for  machine  computation  due  to  the  requirement  of  having  to 
determine  roots  of  polynomials.  In  addition,  the  first  step  would  require  approximating 
all  measured  spectral  matrices  with  rational  spectral  matrices.  This  step  could  also 
be  quite  complicated  and  could  entail  a  serious  loss  with  respect  to  the  amount  of  noise 
power  which  can  be  minimized. 

It  is  difficult  in  general  to  determine  how  well  a  physically  realizable  filter 
performs  relative  to  a  two-sided  filter.  The  two-sided  filter  must  always  be  better 
than  the  physically  realizable  filter  since 

E{So2}  ,  E{§f}  . 

This  follows  from  the  fact  that  the  optimum  weighting  functions  in  the  two-sided  case 
are  subject  to  fewer  constraints  than  those  for  the  physically  realizable  case.  How¬ 
ever,  we  can  obtain  some  results  along  these  lines  by  using  the  theory  of  Toeplitz 

*  12 
forms. 

We  have  already  shown,  cf.  Eq.  (57)  , 

-2  r 

limit  E  {  S  ,  }  =  J 

V  00  J  —  TT 


K 

S 

U,k=l  J  . 


-1 


dx 

2tt 


(90) 
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for  a  fixed  j*.  The  quantity  on  the  right  in  Eq.  (90)  represents  a  lower  bound  for  the 

variance  of  the  processed  noise  and  can  only  be  attained  asymptotically  as  the  memory 

of  the  filter,  (2v  +  1),  increases.  We  obtain  from  Eqs.  (7),  (22),  in  conjunction  with 

12 

some  results  from  theory  of  Toeplitz  forms, 


limit 

v  00 


1 

2v+  1 


U»k=1 


-I 


dx 
2-  • 


(91) 


This  equation  tells  us  what  the  average  processed  noise  variance  should  be,  where 
the  averaging  is  done  with  respect  to  the  time,  j',  at  which  the  signal  is  to  be 
estimated. 
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VI. 


INTERPRETATION  OF  OPTIMUM  MULTIDIMENSIONAL  FILTERING  IN  TERMS 
OF  FREQUENCY  WAVE-NUMBER  SPACE 


We  now  wish  to  show  the  manner  in  which  maximum -likelihood  filtering  can  be 

interpreted  in  the  frequency  wave -number  space.  We  introduce  the  wave -numbers 

k  ,  k  and  the  frequency  wave -number  spectrum  as  follows 
x  y 


K 


f(x\k  ,k  )  =2  E  P.-Wexptil-tx’-k  (x -O-k  (y -y  )] 

y  j,k=l  t=-<*>  JK  x  j  k  y  j  x 


K 


=  Tj  f..{x')exp  (i(k  x.+k  y.)]  exp  [i(k  x  +k  y  ))  , 
j  k~l  Jk  x  j  y  j  x  k  y  k 


(92) 


where  (x.fy.)  and  (x.  ,y, )  are  the  spatial  coordinates  of  the  jtl1  and  kth  seismometers, 
j  j  k  k 

respectively,  measured  with  respect  to  some  reference.  Since  {f  (x’)}  is  a  non- 

JK 

negative -definite  matrix  for  all  x’,  it  follows  from  Eq.  (92)  that  f  is  nonnegative  for 

all  x* ,  k  ,  k  . 
x  y 

We  have 


n  n  _  <}V  dky 

fjk<x')  =  1  J-  f<x',kx,  k  )  exp  l  i(kxx.  +  y.)  ]  exp  [  -i(kx.  +  kyy.)  ]  —  -£■ 
J  —  tt  tt  *'  ’  J 


(93) 


so  that 


^ 

2  tt  2tt  2tt  * 


where 

B(x’,k  ,k  ) 
x  y 


K  ® 

E  E  Mk  l  j*)  exP  l  ^  +  k  x.  +  k  y.) )  • 

j=i  k=-»  j  3  y  3 


(94) 


(95) 
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We  now  consider  only  the  situation  where  j*  =  0,  in  which  case  the  constraint 


Eq.  (1.1)  becomes 


B(x’,0,0)  =  1  . 


(96) 


Thus,  the  maximum -likelihood  filter,  or  minimum -variance  unbiased  estimator,  has 
those  filter  weights  9^ ( k j 0)  such  that  the  triple  integral  over  frequency  wave-number 
space,  in  Eq.  (94)  is  minimized  subject  to  the  constraints  in  Eqs.  (95),  (96).  The 
solution  to  this  problem  would  have  to  be  obtained  in  the  usual  way,  i.e. ,  by  using  the 
recursive  method  of  solution,  assuming,  of  course,  that  v  is  finite,  the  approximate 
frequency  domain  approach  for  two-sided  filters,  if  v  is  large,  or  the  spectral 
matrix  factorization  technique  for  physically  realizable  filters,  when  v  is  1  rge. 

However,  we  can  visualize  to  a  great  extent  what  the  minimizing  solution  for  B 

should  be  in  the  frequency  wave -number  space.  It  is  easily  seen  that  if  K  is  reasonably 

large  and  if  the  spatial  coordinates  <x^ *  y j )  are  reasonably  different  from  (x^,y^),  for 

all  j,k,  the  minimizing  solution  for  B  will  be  such  that  j B |  is  small  whenever  f  is 

large,  except  at  the  origin  of  the  wave -number  plane,  where,  fo~  all  frequencies,  B  is 

equal  to  unity.  Thus,  |b|  would  tend  to  be  needle-shaped,  with  amplitude  unity,  at 

the  origin  of  the  k^-  k^  plane,  for  all  values  of  the  frequc  icy  variable  x’,  and  would 

tend  to  have  sidelobes  in  those  parts  of  the  k  -k  plane  where  f  would  be  small,  for  a 

x  y 

given  x\  In  general,  the  positions  and  amplitudes  of  these  lobes  would  be  frequency 
dependent. 


35 


vn.  CONCLUSIONS 


It  has  been  shown  that,  basically,  the  reason  for  the  equivalence  of  the 
maximum  -likelihood  and  minimum  -variance  unbiased  estimators  is  that  they  are 
based  on  a  conditional  expectation.  The  martingale  property  of  conditional  expectation 
then  assures  that  the  asymptotic  properties  of  these  estimators  are  well  defined. 

An  asymptotically  optimum  frequency  domain  synthesis  procedure  has  been 
presented  which  requires  much  less  computing  time  than  the  exact  recursive  procedure. 
However,  the  frequency  domain  synthesis  technique  is  valid  only  for  two-sided  filters, 
while  the  recursive  procedure  is  always  valid.  If  two-sided  filters  are  to  be  used  due 
to  their  inherently  better  capability  to  suppress  the  noise,  then  the  frequency  domain 
approach  is  superior,  assuming  that  v  is  large  enough  for  the  asymptotic  results  to 
hold. 

Synthesis  procedures  for  physically  realizable  filters,  in  the  asymptotic  case, 
have  also  been  presented.  These  procedures  were  based  on  a  spectral  matrix  factori¬ 
zation  technique,  and  are  not  well  adapted  to  machine  computation. 
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